First-Order Methods

19.8 μs

Spoiler alert: Some methods are trash!

13.8 μs

68.5 μs

37 μs

169

139 μs

0.27

72 μs

0.82

82.8 μs

16.1 ms

67.6 ms

Gradient Descent

9.8 μs

xxxxxxxxxx
 
using Calculus

231 ms

We can choose the direction of steepest descent as the descent direction d, and this is the opposite of the gradient.

11.6 μs

$d^{(k)} = - \frac{g^{(k)}}{| | g^{(k)} | |}$

20.1 μs

It is proven (proof left to the reader as excercise) that each direction is orthogonal to the previous.

12.2 μs

xxxxxxxxxx
 
using Base.MathConstants

4.9 μs

fibonacci_search (generic function with 1 method)

168 μs

bracket_minimum (generic function with 2 methods)

 
function bracket_minimum(f::Function, x=0.0; s=1e-12, k=2.0)
    try
        a, ya = x, f(x)
        global b, yb = x+s, f(x+s)
        if yb > ya #If we are not going downhill
            c, yc = x-s, f(x-s)
            s = -s
        else
            c, yc = x+k*s, f(x+k*s) 
        end
        if ya > yb < yc # Definition of unimodal
            return (a, c)
        end
        s = s*k
    catch e
        println(e)
        return (a, c)
    end
    return bracket_minimum(f, b; s=s, k=k)
end

120 μs

minimize (generic function with 1 method)

x
 
minimize(f::Function, a::Number, b::Number)::Number = fibonacci_search(f, a, b, 5)

34.1 μs

optimize_line_search (generic function with 1 method)

xxxxxxxxxx
 
function optimize_line_search(f, x, d)::Float64
    objective(α::Number)::Number = f(x + α*d)
    a, b = bracket_minimum(objective; k=1.05)
    α = minimize(objective, a, b)
    α
end

71.5 μs

gradient_descent_mixed (generic function with 1 method)

xxxxxxxxxx
 
function gradient_descent_mixed(f, ∇f, x, k)::Array{Tuple, 1}
    points::Array{Tuple} = Array{Tuple}(undef, k+1)
    points[1] = tuple(x...)
    for i in 1:k
        d = -∇f(x)
        d = d/sqrt(sum(2 .^ d))
        α = optimize_line_search(f, x, d)
        x = x + α*d
        points[i+1] = tuple(x...)
    end
    points
end

49.4 μs

f₁ (generic function with 2 methods)

xxxxxxxxxx
 
begin
    f₁(x, y)::Number = (1-x)^2 + 5(y - x^2)^2
    f₁(x̄)::Number = f₁(x̄...)
end

60.9 μs

xxxxxxxxxx
 
using Plots

1.1 ms

xxxxxxxxxx
 
using PlutoUI

184 ms

45.4 ms

-0.96

36.5 ms

1.69

78.6 μs

x
 
begin
    x₁ = range(-3, stop=3, length=200)
    y₁ = x₁
    color₁ = cgrad([:grey, :blue])
    ∇f₁(x̄) = Calculus.gradient(f₁, x̄)
    points₁ = [x for x in gradient_descent_mixed(f₁, ∇f₁, [x₁ᵢ, y₁ᵢ], n₁)]
    plot₁ = contour(x₁, y₁, f₁, c=color₁, levels=700 .^ (range(-.2,stop=1,length=14)))
    scatter!((1.0, 1.0), markersize=3, markerstrokewidth=0, c=:cyan, label="Optimal point")
    scatter!((x₁ᵢ, y₁ᵢ), markersize=2, c=:black, label="Initial point")
    plot!(points₁, c=:red, label="Descent")
    plot₁
end

138 ms

The optimal point is near (0.9605721406277097, 0.9196323117617097)

14.7 μs

The optimized value is 0.00160157399060238

14.7 μs

xxxxxxxxxx
 
abstract type DescentMethod end

10.2 μs

xxxxxxxxxx
 
begin
    struct GradientDescent <: DescentMethod
        α
    end
end

1.6 ms

init! (generic function with 1 method)

xxxxxxxxxx
 
init!(M::GradientDescent, f, ∇f, x) = M

31.7 μs

step! (generic function with 1 method)

xxxxxxxxxx
 
function step!(M::GradientDescent, f, ∇f, x)::Array{Float64, 1}
    α, g = M.α, ∇f(x)
    return x - α*g
end

34.2 μs

optimize! (generic function with 1 method)

xxxxxxxxxx
 
function optimize!(M::DescentMethod, f, ∇f, x; k=30)::Array{Tuple}
    init!(M, f, ∇f, x)
    points::Array{Tuple} = Array{Tuple}(undef, k+1)
    for i in 1:k
        points[i] = tuple(x...)
        x = step!(M, f, ∇f, x)
    end
    points[end] = tuple(x...)
    points
end

72 μs

f₂

f₁ (generic function with 2 methods)

xxxxxxxxxx
 
f₂ = f₁

1.5 μs

95.5 μs

-0.96

83.3 μs

1.69

90.9 μs

37 ms

The optimal point is near (0.7414365199553092, 0.6976297924702967) 👀

9.4 μs

The optimized value is 0.17622960698092363 😒🤢

12.6 μs

Gradient descent (alone) is not miraculous 😢

9.9 μs

Conjugate Gradient

10.5 μs

Gradient descent is 🗑️ in narrow valleys. This method uses a quadratic approximation of the local function to find optimal points.

9.7 μs

Conjugate Gradient Descent is GUARANTEED to optimize a n-dimensional quadratic function in n steps.

10.4 μs

This method can be applied to nonquadratic functions. Smooth functions behave quadratic close to a local minimum (oh-oh).

9.7 μs

This means that we need to get close to the minimum in order to obtain benefits. 🗑️

8.4 μs

This method uses information from last gradients to improve:

9.8 μs

$d^{(k + 1)} = - g^{(k + 1)} + β^{(k)} d^{(k)}$

12.6 μs

Fletcher-Reeves: $β^{(k)} = \frac{(g^{(k)})^{2}}{(g^{(k - 1)})^{2}}$

21.1 μs

Polak-Ribiere: $β^{(k)} = \frac{g^{(k) T} (g^{(k)} - g^{(k - 1)})}{(g^{(k - 1)})^{2}}$

15.6 μs

👀 For guaranteed convergence on the Polak-Ribiere method we need to reset β to 0 if it goes negative.

11.2 μs

ConjugateGradientDescent

xxxxxxxxxx
 
begin
    mutable struct ConjugateGradientDescent <: DescentMethod
        d::Array{Number, 1}
        g::Array{Number, 1}
    end
    function ConjugateGradientDescent()
        return ConjugateGradientDescent(Number[], Number[])
    end
end

4.2 ms

init! (generic function with 2 methods)

xxxxxxxxxx
 
function init!(M::ConjugateGradientDescent, f, ∇f, x)
    M.g = ∇f(x)
    M.d = -M.g
    return M
end

34.4 μs

xxxxxxxxxx
 
using LinearAlgebra

15.7 ms

step! (generic function with 10 methods)

xxxxxxxxxx
 
function step!(M::ConjugateGradientDescent, f, ∇f, x)::Array{Float64, 1}
    d, g = M.d, M.g
    g′ = ∇f(x)
    β = (g′ ⋅ (g′ - g))/(g ⋅ g)
    β = β*(β>0)
    d′ = -g′ + β*d
    α = optimize_line_search(f, x, d′)
    x′ = x + α*d′
    M.d, M.g = d′, g′
    return x′
end

55.8 μs

f₃

f₁ (generic function with 2 methods)

xxxxxxxxxx
 
f₃ = f₁

1.9 μs

57.8 μs

-0.96

80.3 μs

1.69

93 μs

34.3 ms

The optimal point is near (0.9926081864948202, 0.9847483130588324) 🌟

17.9 μs

The optimized value is 5.6004977270041056e-5 🥳

20.3 μs

Momentum

8.7 μs

The idea is to simulate a force (that updates velocity) to avoid slow movement in a nearly flat surface (gradients with small magnitude will fail miserably).

7.4 μs

$v^{(k + 1)} = β v^{(k)} - α g^{(k)}$

10.1 μs

$x^{(k + 1)} = x^{(k)} + v^{(k + 1)}$

8.8 μs

Momentum

xxxxxxxxxx
 
begin
    mutable struct Momentum <: DescentMethod
        α::Number
        β::Number
        v::Array{Number, 1}
    end
    Momentum(α, β) = Momentum(α, β, Number[])
end

3.1 ms

init! (generic function with 3 methods)

xxxxxxxxxx
 
function init!(M::Momentum, f, ∇f, x)
    M.v = zeros(length(x))
end

32.1 μs

step! (generic function with 3 methods)

xxxxxxxxxx
 
function step!(M::Momentum, f, ∇f, x)::Array{Float64, 1}
    α, β, v, g = M.α, M.β, M.v, ∇f(x)
    v[:] = β*v - α*g #Update velocity
    return x + v
end

52.1 μs

f₄ (generic function with 2 methods)

xxxxxxxxxx
 
begin
    f₄(x, y)::Number = (1-x)^2 + 100(y - x^2)^2
    f₄(x̄)::Number = f₄(x̄...)
end

58.1 μs

77.2 μs

-0.58

67.5 μs

0.36

6.3 ms

xxxxxxxxxx
 
begin
    x₄ = range(-1.5, stop=1.25, length=200)
    y₄ = range(-.45, stop=1.5, length=200)
    color₄ = cgrad([:yellow,:green, :blue, :red, :purple])
    ∇f₄(x̄) = Calculus.gradient(f₄, x̄)
    points₄ = [x for x in optimize!(Momentum(.005, 1), f₄, ∇f₄, [x₄ᵢ, y₄ᵢ]; k=n₄)]
    contour(x₄, y₄, f₄, c=color₄, levels=2500 .^ (range(-.5,stop=1,length=25)), legend=false)
    scatter!((1.0, 1.0), markersize=3, markerstrokewidth=0, c=:cyan, label="Optimal point")
    scatter!((x₄ᵢ, y₄ᵢ), markersize=4, c=:black, label="Initial point")
    plot!(points₄, c=:red, label="Descent")
    scatter!(points₄, c=:green, markersize=2.5, markerstrokewidth=0, label="Path")
end

9.1 ms

The optimal point is near (1.0263932677227146, 1.304148459285668) 🤷‍♀️

13.2 μs

The optimized value is 6.28400683251055 🤦‍♂️

10.5 μs

It's like the point is falling straight to the minimum 😲

7.8 μs

But... it might go too fast

7.5 μs

Nesterov Momentum

6.7 μs

If it goes too fast try predicting how fast it will go the next iteration...

10 μs

$v^{(k + 1)} = β v^{(k)} - α \nabla f (x^{(k)} + β v^{(k)})$

9.3 μs

step! (generic function with 4 methods)

xxxxxxxxxx
 
begin
    mutable struct NesterovMomentum <: DescentMethod
        α::Number
        β::Number
        v::Array{Number, 1}
    end
    function NesterovMomentum(α::Number, β::Number)
        NesterovMomentum(α, β, Number[])
    end
    function init!(M::NesterovMomentum, f, ∇f, x)
        M.v = zeros(length(x))
        return M
    end
    function step!(M::NesterovMomentum, f, ∇f, x)::Array{Float64, 1}
        α, β, v = M.α, M.β, M.v
        next_g = β*v
        v[:] = next_g - α*∇f(x + next_g)
        return x + v
    end
end

2.3 ms

f₅

f₄ (generic function with 2 methods)

xxxxxxxxxx
 
f₅ = f₄

3 μs

1000

xxxxxxxxxx
 
@bind n₅ Slider(1:1000, default=42, show_value=true)

68.9 μs

-0.77

xxxxxxxxxx
 
@bind x₅ᵢ Slider(-1.5:.01:1.5, default=-0.77, show_value=true)

64.8 μs

1.42

xxxxxxxxxx
 
@bind y₅ᵢ Slider(-1.5:.01:1.5, default=1.42, show_value=true)

110 μs

x
 
begin
    x₅ = range(-1.5, stop=1.25, length=200)
    y₅ = range(-.45, stop=1.5, length=200)
    color₅ = cgrad([:yellow,:green, :blue, :red, :purple])
    ∇f₅(x̄) = Calculus.gradient(f₅, x̄)
    points₅ = [x for x in optimize!(NesterovMomentum(.0009, .95), f₅, ∇f₅, [x₅ᵢ, y₅ᵢ]; k=n₅)]
    contour(x₅, y₅, f₅, c=color₅, levels=2500 .^ (range(-.5,stop=1,length=25)), legend=false)
    scatter!((1.0, 1.0), markersize=3, markerstrokewidth=0, c=:cyan, label="Optimal point")
    scatter!((x₅ᵢ, y₅ᵢ), markersize=4, c=:black, label="Initial point")
    plot!(points₅, c=:red, label="Descent")
    scatter!(points₅, c=:green, markersize=1.5, markerstrokewidth=0, label="Path")
end

15.6 ms

The optimal point is near (0.9999678362763855, 0.9999355449078215) 🤔

12.2 μs

The optimized value is 1.03616095695321e-9 🤔🤔

17.7 μs

It's still going too fast...

7.4 μs

Adagrad

7.5 μs

Adaptative subgradient method

8.7 μs

This method updates learning rate for each component.

8.3 μs

$x_{i}^{(k + 1)} = x_{i}^{(k)} - \frac{α}{ϵ + \sqrt{s_{i}^{(k)}}} g_{i}^{(k)}$

9 μs

$s_{i}^{(k)} = \sum_{j = 1}^{k} (g_{i}^{(j)})^{2}$

7.8 μs

ϵ is a small value (1e-8) to prevent division by zero.

7.4 μs

α is usually set to .01, because of all the operations it doesn't matter too much

8.9 μs

Problem: As we can see, the components of S grow over time, making the learning rate decrease (and become infinitesimally small even before convergence).

8 μs

step! (generic function with 5 methods)

xxxxxxxxxx
 
begin
    mutable struct Adagrad <: DescentMethod
        α::Number
        ϵ::Number
        s::Array{Float64, 1}
    end
    Adagrad(α::Number, ϵ::Number) = Adagrad(α, ϵ, Float64[])
    Adagrad(α::Number) = Adagrad(α, 1e-8)
    function init!(M::Adagrad, f, ∇f, x)
        M.s = zeros(length(x))
        return M
    end
    function step!(M::Adagrad, f, ∇f, x)::Array{Float64, 1}
        α, ϵ, s, g = M.α, M.ϵ, M.s, ∇f(x)
        s[:] += g .* g
        return x - α*g ./ (sqrt.(s) .+ ϵ)
    end
end

3.5 ms

f₆

f₄ (generic function with 2 methods)

xxxxxxxxxx
 
f₆ = f₄

2.3 μs

107 μs

-0.77

112 μs

14.6 ms

68.3 ms

The optimal point is near (-0.3043008040570129, 0.09737789529620958) 🤣

14.9 μs

The optimized value is 1.7034843912261084 🐌🐌🐌

14.6 μs

Really slow convergence...

7.7 μs

RMSProp

8 μs

It comes to the rescue of Adagrad, solving the problem of decreasing learning rate.

11.3 μs

⊙ = \odot (element wise multiplication)

8.2 μs

$s ̂^{(k + 1)} = γ s ̂^{(k)} + (1 - γ) (g^{(k)} ⊙ g^{(k)})$

8.3 μs

γ is tipically close to 0.9

8.1 μs

$R M S (g_{i}) = s ̂_{i}^{(k + 1)}$

8.4 μs

step! (generic function with 6 methods)

xxxxxxxxxx
 
begin
    mutable struct RMSProp <: DescentMethod
        α::Number
        γ::Number
        ϵ::Number
        ŝ::Array{Float64, 1}
    end
    RMSProp(α::Number, γ::Number, ϵ::Number) = RMSProp(α, γ, ϵ, Float64[])
    RMSProp(α::Number, γ::Number) = RMSProp(α, γ, 1e-8)
    RMSProp(α::Number) = RMSProp(α, .9)
    function init!(M::RMSProp, f, ∇f, x)
        M.ŝ = zeros(length(x))
        return M
    end
    function step!(M::RMSProp, f, ∇f, x)
        α, γ, ϵ, ŝ, g = M.α, M.γ, M.ϵ, M.ŝ, ∇f(x)
        ŝ[:] = γ*ŝ + (1-γ)*(g .* g)
        return x - α*g ./ (sqrt.(ŝ) .+ ϵ)
    end
end

4 ms

f₇

f₄ (generic function with 2 methods)

xxxxxxxxxx
 
f₇ = f₄

4.4 μs

294

124 μs

-1.36

141 μs

0.13

146 μs

x
 
begin
    x₇ = range(-1.5, stop=1.5, length=200)
    y₇ = range(-.45, stop=1.5, length=200)
    color₇ = cgrad([:yellow,:green, :blue, :red, :purple])
    ∇f₇(x̄) = Calculus.gradient(f₇, x̄)
    points₇ = [x for x in optimize!(RMSProp(.04, .999, 1), f₇, ∇f₇, [x₇ᵢ, y₇ᵢ]; k=n₇)]
    contour(x₇, y₇, f₇, c=color₇, levels=2500 .^ (range(-.5,stop=1,length=25)), legend=false)
    scatter!((1.0, 1.0), markersize=3, markerstrokewidth=0, c=:cyan, label="Optimal point")
    scatter!((x₇ᵢ, y₇ᵢ), markersize=4, c=:black, label="Initial point")
    plot!(points₇, c=:red, label="Descent")
    scatter!(points₇, c=:green, markersize=1, markerstrokewidth=0, label="Path")
end

37.1 ms

Adadelta

8.6 μs

Removes the learning rate parameter entirely (the same authors of RMSProp did this) 😲😲

7.1 μs

$x_{i}^{(k + 1)} = x_{i}^{(k)} - \frac{R M S (Δ x_{i})}{ϵ + R M S (g_{i})}$

6.7 μs

step! (generic function with 7 methods)

xxxxxxxxxx
 
begin
    mutable struct Adadelta <: DescentMethod
        γs::Number
        γx::Number
        ϵ::Number
        ŝ::Array{Float64, 1}
        u::Array{Float64, 1}
    end
    Adadelta(γs::Number, γx::Number, ϵ::Number) = Adadelta(γs, γx, ϵ, Float64[], Float64[])
    Adadelta(γs::Number, γx::Number) = Adadelta(γs, γx, 1e-2)
    function init!(M::Adadelta, f, ∇f, x)
        M.ŝ = zeros(length(x))
        M.u = zeros(length(x))
        return M
    end
    function step!(M::Adadelta, f, ∇f, x)::Array{Float64, 1}
        γs, γx, ϵ, ŝ, u, g = M.γs, M.γx, M.ϵ, M.ŝ, M.u, ∇f(x)
        ŝ[:] = γs*ŝ + (1-γs)*(g .* g)
        Δx = -(sqrt.(u) .+ ϵ) ./ (sqrt.(ŝ) .+ ϵ) .* g
        u[:] = γx*u + (1-γx)*(Δx .* Δx)
        return x + Δx
    end
end

5.1 ms

f₈

f₄ (generic function with 2 methods)

xxxxxxxxxx
 
f₈ = f₄

2.6 μs

577

x
 
@bind n₈ Slider(1:1000, default=102, show_value=true)

67.6 μs

-1.06

110 μs

0.87

105 μs

xxxxxxxxxx
 
begin
    x₈ = range(-1.5, stop=1.5, length=200)
    y₈ = range(-.45, stop=1.5, length=200)
    color₈ = cgrad([:yellow,:green, :blue, :red, :purple])
    ∇f₈(x̄) = Calculus.gradient(f₈, x̄)
    points₈ = [x for x in optimize!(Adadelta(.9, .9999, 3e-2), f₈, ∇f₈, [x₈ᵢ, y₈ᵢ]; k=n₈)]
    contour(x₈, y₈, f₈, c=color₈, levels=2500 .^ (range(-.5,stop=1,length=25)), legend=false)
    scatter!((1.0, 1.0), markersize=3, markerstrokewidth=0, c=:cyan, label="Optimal point")
    scatter!((x₈ᵢ, y₈ᵢ), markersize=4, c=:black, label="Initial point")
    plot!(points₈, c=:red, label="Descent")
    scatter!(points₈, c=:green, markersize=1, markerstrokewidth=0, label="Path")
end

45.7 ms

Adam

4.6 μs

Adaptative moment estimation method

4.9 μs

Stores exponentially decaying squared gradient (RMSProp and Adadelta), but also decaying gradient like momentum.

4.6 μs

Initializing the gradient and squared gradient to zero introduces a bias. (Good defaults are α=.001, γᵥ=.9, γₛ=.999 and ϵ=1e-8)

5 μs

The equations for Adam are:

4.7 μs

Biased decaying momentum: $v^{(k + 1)} = γ_{v} v^{(k)} + (1 - γ_{v}) g^{(k)}$

4 ms

Biased decaying sq. gradient: $s^{(k + 1)} = γ_{s} s^{(k)} + (1 - γ_{s}) (g^{(k)} ⊙ g^{(k)})$

10.9 μs

Corrected decaying momentum: $v ̂^{(k + 1)} = v^{(k + 1)} / (1 - γ_{v}^{k})$

8.5 μs

Corrected decaying sq. gradient: $s ̂^{(k + 1)} = s^{(k + 1)} / (1 - γ_{s}^{k})$

10.9 μs

Next iterate: $x^{(k + 1)} = x^{(k)} - α v ̂^{(k + 1)} / (ϵ + \sqrt{s ̂^{(k + 1)}})$

12 μs

step! (generic function with 8 methods)

xxxxxxxxxx
 
begin
    mutable struct Adam <: DescentMethod
        α::Float64
        γᵥ::Float64
        γₛ::Float64
        ϵ::Float64
        k::Integer
        v::Array{Float64, 1}
        s::Array{Float64, 1}
    end
    Adam(α::Number, γᵥ::Number, γₛ::Number, ϵ::Number) = Adam(α, γᵥ, γₛ, ϵ, 0, Float64[], Float64[])
    Adam(α::Number, γᵥ::Number, γₛ::Number) = Adam(α, γᵥ, γₛ, 1e-8)
    Adam(α::Number, γᵥ::Number) = Adam(α, γᵥ, .999)
    Adam(α::Number) = Adam(α, .9)
    
    function init!(M::Adam, f, ∇f, x)
        M.k = 0
        M.v = zeros(length(x))
        M.s = zeros(length(x))
        return M
    end
    
    function step!(M::Adam, f, ∇f, x)
        α, γᵥ, γₛ, ϵ, k = M.α, M.γᵥ, M.γₛ, M.ϵ, M.k
        s, v, g = M.s, M.v, ∇f(x)
        v[:] = γᵥ*v + (1-γᵥ)*g
        s[:] = γₛ*s + (1-γₛ)*(g .* g)
        M.k = k += 1 #Updates k at the same time
        v̂ = v/(1-γᵥ^k)
        ŝ = s/(1-γₛ^k)
        return x - (α*v̂)./(sqrt.(ŝ) .+ ϵ)
    end
end

6.3 ms

f₉

f₄ (generic function with 2 methods)

xxxxxxxxxx
 
f₉ = f₄

4.2 μs

285

x
 
@bind n₉ Slider(1:500, default=222, show_value=true)

61.8 μs

-0.98

105 μs

0.84

107 μs

x
 
begin
    x₉ = range(-1.5, stop=1.5, length=200)
    y₉ = range(-.45, stop=1.5, length=200)
    color₉ = cgrad([:yellow,:green, :blue, :red, :purple])
    ∇f₉(x̄) = Calculus.gradient(f₉, x̄)
    points₉ = [x for x in optimize!(Adam(.3, .95, .99), f₉, ∇f₉, [x₉ᵢ, y₉ᵢ]; k=n₉)]
    contour(x₉, y₉, f₉, c=color₉, levels=2500 .^ (range(-.5,stop=1,length=25)), legend=false)
    scatter!((1.0, 1.0), markersize=3, markerstrokewidth=0, c=:cyan, label="Optimal point")
    scatter!((x₉ᵢ, y₉ᵢ), markersize=4, c=:black, label="Initial point")
    plot!(points₉, c=:red, label="Descent")
    scatter!(points₉, c=:green, markersize=1.5, markerstrokewidth=0, label="Path")
end

50.9 ms

Hypergradient Descent

7.1 μs

This methods are too sensitive to the learning rate... let's optimize it first

8.3 μs

We applied gradient descent to the learning reate 🤯🤯

8.5 μs

step! (generic function with 9 methods)

xxxxxxxxxx
 
begin
    mutable struct HyperGradientDescent <: DescentMethod
        α₀::Float64 #Initial learning rate
        μ::Float64 #Inception learning rate
        α::Float64 #Current learning rate
        g_prev::Array{Float64, 1}
    end
    
    HyperGradientDescent(α₀::Number, μ::Number) = HyperGradientDescent(α₀, μ, α₀, Float64[])
    
    function init!(M::HyperGradientDescent, f, ∇f, x)
        M.α = M.α₀
        M.g_prev = zeros(length(x))
        return M
    end
    
    function step!(M::HyperGradientDescent, f, ∇f, x)::Array{Float64, 1}
        α, μ, g, g_prev = M.α, M.μ, ∇f(x), M.g_prev
        α = α + μ*(g ⋅ g_prev)
        M.g_prev, M.α = g, α
        return x - α*g
    end
end

3 ms

f₁₀

f₄ (generic function with 2 methods)

xxxxxxxxxx
 
f₁₀ = f₄

3.1 μs

102 μs

-0.77

128 μs

101 μs

xxxxxxxxxx
 
begin
    x₁₀ = range(-1.5, stop=1.5, length=200)
    y₁₀ = range(-.45, stop=1.5, length=200)
    color₁₀ = cgrad([:yellow,:green, :blue, :red, :purple])
    ∇f₁₀(x̄) = Calculus.gradient(f₁₀, x̄)
    points₁₀ = [x for x in optimize!(HyperGradientDescent(.005, 1e-8), f₁₀, ∇f₁₀, [x₁₀ᵢ, y₁₀ᵢ]; k=n₁₀)]
    contour(x₁₀, y₁₀, f₁₀, c=color₁₀, levels=2500 .^ (range(-.5,stop=1,length=25)), legend=false)
    scatter!((1.0, 1.0), markersize=3, markerstrokewidth=0, c=:cyan, label="Optimal point")
    scatter!((x₁₀ᵢ, y₁₀ᵢ), markersize=4, c=:black, label="Initial point")
    plot!(points₁₀, c=:red, label="Descent")
    scatter!(points₁₀, c=:green, markersize=1.5, markerstrokewidth=0, label="Path")
end

39.4 ms

Let's do the same to Nesterov momentum

6.6 μs

step! (generic function with 10 methods)

xxxxxxxxxx
 
begin
    mutable struct HyperNesterovMomentum <: DescentMethod
        α₀::Float64 #Initial learning rate
        μ::Float64 #Inception learning rate
        β::Float64
        v::Array{Float64, 1}
        α::Float64
        g_prev::Array{Float64, 1}
    end
    
    HyperNesterovMomentum(α₀::Number, μ::Number, β::Number) = HyperNesterovMomentum(α₀, μ, β, Float64[], α₀, Float64[])
    
    function init!(M::HyperNesterovMomentum, f, ∇f, x)
        M.α = M.α₀
        M.v = zeros(length(x))
        M.g_prev = zeros(length(x))
        return M
    end
    
    function step!(M::HyperNesterovMomentum, f, ∇f, x)::Array{Float64, 1}
        α, β, v = M.α, M.β, M.v
        g, μ, g_prev = ∇f(x), M.μ, M.g_prev
        next_g = β*v
        α = α + μ*(g ⋅ (-g_prev - β*v))
        v[:] = β*v + g
        M.α, M.g_prev = α, g_prev
        return x - α*(g + β*v)
    end
end

2.5 ms

f₁₁

f₄ (generic function with 2 methods)

xxxxxxxxxx
 
f₁₁ = f₄

1.8 μs

200

77.3 μs

-0.77

76.3 μs

64.9 μs

xxxxxxxxxx
 
begin
    x₁₁ = range(-1.5, stop=1.5, length=200)
    y₁₁ = range(-.45, stop=1.5, length=200)
    color₁₁ = cgrad([:yellow,:green, :blue, :red, :purple])
    ∇f₁₁(x̄) = Calculus.gradient(f₁₁, x̄)
    points₁₁ = [x for x in optimize!(HyperNesterovMomentum(.5e-3, 3e-9, .97), f₁₁, ∇f₁₁, [x₁₁ᵢ, y₁₁ᵢ]; k=n₁₁)]
    contour(x₁₁, y₁₁, f₁₁, c=color₁₀, levels=2500 .^ (range(-.5,stop=1,length=25)), legend=false)
    scatter!((1.0, 1.0), markersize=3, markerstrokewidth=0, c=:cyan, label="Optimal point")
    scatter!((x₁₁ᵢ, y₁₁ᵢ), markersize=4, c=:black, label="Initial point")
    plot!(points₁₁, c=:red, label="Descent")
    scatter!(points₁₁, c=:green, markersize=1.5, markerstrokewidth=0, label="Path")
end

19 ms

plot_method! (generic function with 1 method)

x
 
function plot_method!(M::DescentMethod, f, ∇f, x, k; curve=:red, path=:green)
    methodname = string(nameof(typeof(M)))
    points = [x for x in optimize!(M, f, ∇f, x; k=k)]
    plot!(points, c=curve, label="$(methodname)")
    scatter!(points, c=path, markersize=1, markerstrokewidht=0, label=nothing)
end

125 μs

Exercises

10 μs

Exercise 5.1. Compute the gradient of $x^{T} A x + b^{T} x$ when A is symmetric

20 μs

Answer: By decomposing the matrix as a sum (really long) and then replacing the sum as a matrix form we get: $\nabla f = 2 A x + b$

16.8 μs

Exercise 5.2. Apply gradient descent with unit step to f(x)=x⁴, compute two iterations.

10.3 μs

∇fₑ (generic function with 2 methods)

x
 
begin
    fₑ(x::Float64)::Array{Float64, 1} = [x^4]
    fₑ(x̄::Array{Float64, 1})::Array{Float64, 1} = fₑ(x̄[1])
    ∇fₑ(x::Float64)::Array{Float64, 1} = [4*x^3]
    ∇fₑ(x̄::Array{Float64, 1})::Array{Float64, 1} = ∇fₑ(x̄[1])
end

136 μs

x
 
begin
    xₑ = range(-.01, .5, length=50)
    plot(xₑ, (x) -> fₑ(x)[1], label="Function")
    pointsₑ = [x[1] for x in optimize!(GradientDescent(1), fₑ, ∇fₑ, [.4]; k=2)]
    plot!(pointsₑ, (x) -> fₑ(x)[1], c=:red, label="Gradient descent")
    scatter!(pointsₑ, (x) -> fₑ(x)[1], c=:black, markersize=1, markerstrokewidht=0, label=nothing)
end

623 ms

Exercise 5.3. Apply one step of gradient descent to f(x)=e^x+e^-x from x =10 with both a unit step and with exact line search.

11.5 μs

∇fₘ (generic function with 2 methods)

x
 
begin
    fₘ(x::Float64)::Array{Float64, 1} = [exp(x)+exp(-x)]
    fₘ(x̄::Array{Float64, 1})::Array{Float64, 1} = fₘ(x̄[1])
    ∇fₘ(x::Float64)::Array{Float64, 1} = Calculus.derivative(fₘ, x)
    ∇fₘ(x̄::Array{Float64, 1})::Array{Float64, 1} = ∇fₘ(x̄[1])
end

55.5 μs

Answer: If someone can do this without the gradient exploting to infinity really quick please help.

25 μs

Exercise 5.4. The conjugate gradient method can be used to find a search direction d when a local quadratic model of a function is available at the current point. With d as search direction, let the model be $q (d) = d^{T} H d + b^{T} d + c$ for a symmetric matrix H, What is the Hessian in this case? What is the gradient of q when d=0? What can go wrong if the conjugate gradient method is applied to the quadratic model to get the search direction d?

12.4 μs

Answer: We knew that: $\nabla q = 2 H d$ , then the gradient (gradient of gradient) will be: $\nabla^{2} q = 2 H$

8.4 μs

When d = 0 the gradient will also be 0. And if this happens we will be dividing by infinitesimal values when computing the next search direction. (This happened to me when using too many iterations)

7.4 μs

Exercise 5.5. How is Nesterov momentum an improvement over momentum?

8.2 μs

Answer: The problem with momentum is that it doesn't slow down, and Nesterov momentum uses the calculation of the gradient at the next step (if it overshoots the gradient will make it go back, correcting the step).

8.9 μs

Exercise 5.6. In what way is the conjugate gradient method an improvement over steepest descent?

9.6 μs

Answer: Gradient descent is really slow near valleys (because the next direction will always be orthogonal to the previous, causing a zig-zag movement), that is corrected when the previous gradient contributes to the next direction (making it more like a straight path to the minimum).

8 μs

Exercise 5.7. In the conjugate gradient descent, what is the normalized descent direction at the first iteration for the function $f (x, y) = x^{2} + x y + y^{2} + 5$ when initialized at (x,y) = (1, 1)? What is the resulting point after two steps of the conjugate gradient method?

10.1 μs

∇fₚ (generic function with 1 method)

x
 
begin
    fₚ(x, y) = x^2 + x*y + y^2 + 5
    fₚ(x̄) = fₚ(x̄...)
    ∇fₚ(x̄) = Calculus.gradient(fₚ, x̄)
end

42.9 μs

x
 
begin
    plotly()
    xₚ = range(-3, 3, length=40)
    yₚ = range(-3, 3, length=40)
    pₚ = contour(xₚ, yₚ, fₚ)
    plot_method!(ConjugateGradientDescent(), fₚ, ∇fₚ, [1, 1], 2)
    gr()
    pₚ
end

2.7 ms

x
 
begin
    Mₚ = ConjugateGradientDescent()
    init!(Mₚ, fₚ, ∇fₚ, [1, 1])
    xₚ₀ = [1, 1]
    xₚ₁ = step!(Mₚ, fₚ, ∇fₚ, xₚ₀)
    r₁, r₂ = Mₚ.d[1], Mₚ.d[2]
    xₚ₂ = step!(Mₚ, fₚ, ∇fₚ, xₚ₁)
end;

1.2 ms

Answer: The direction is [-2.9999999999239875, -2.9999999999239875], (Almost the same on each component), And the point is [1.0056882026629292e-6, 1.0056882026629292e-6] after two iterations (expected because this function is polinomial of order two and the conjugate gradient should optimize it in two steps).

14.4 μs

Exercise 5.8. We have a polynomial function f such that f(x) > 2 for all x in three-dimensional Euclidean space. Suppose we are using steepest descent with step lengths optimized at each step, and we want to find a local minimum of f. If our unnormalized descent direction is [1,2,3] at step k, is it possible for our unnormalized descent direction at step k+1 to be [0,0,-3]? Why or why not?

8.5 μs

Answer: We know that in the steepest descent method the descent direction in every iteration is orthogonal to the next direction. $[1, 2, 3] \cdot [0, 0, 3] \neq 0$ Then the unnormalized descent direction at step k+1 cannot be [0,0,-3]. (TODO: Check the restriction for f(x) > 2)

8.1 μs